Importing necessary libraries

In [1]:
import numpy as np
from scipy.io import wavfile
import matplotlib.pyplot as plt 
from scipy.interpolate import interp1d
from scipy.signal import tf2zpk, find_peaks, freqz
from scipy.io.wavfile import write

Please see

Sounds are generated and stored in .wav file in the same folder

Please see : Qn7 : /aa/ sound resynthesized - $aaReconstructed.wav$

Qn7 : Pitch variation - $aaReconstructedPitchVariationLow.wav$, $aaReconstructedPitchVariationHigh.wav$

Bonus : /ss/ sound resynthesized - $ssReconstructed.wav$

Question 1

Applying pre-emphasis to the signal

We can do this by passing the signal througha filter : $$H(z) = (1 - \alpha z^{-1}) \text{ with $\alpha$ very close to 1, typically = 0.95}$$ $$\implies y[n] = x[n] - \alpha x[n-1]$$

In [2]:
pathToaaSound = r"aa.wav" #path to the "pani" sound of machali
sampleRate, data = wavfile.read(pathToaaSound) #storing the used sample rate and the wav file in variables for analysis
data = data/(32767.0)#Normalising the data
data
Out[2]:
array([ 3.56761376e-02,  6.76900540e-02, -1.15970336e-03, -7.36716819e-02,
       -1.66692099e-01, -1.71269875e-01, -1.23935667e-01, -4.41297647e-02,
        3.31431013e-02,  2.82601398e-02, -1.73955504e-03, -5.22782067e-02,
       -5.78936125e-02, -3.38755455e-03,  7.37327189e-02,  1.55125584e-01,
        1.54820399e-01,  1.48045289e-01,  1.21890927e-01,  1.58085879e-01,
        2.27454451e-01,  3.13638722e-01,  3.98754845e-01,  2.09662160e-01,
       -9.47599719e-02, -2.49916074e-01, -3.69609668e-01, -2.82631916e-01,
        1.22684408e-02,  1.87810907e-01,  2.54982147e-01,  1.93914609e-01,
        1.08798486e-01, -3.26242866e-02, -4.34278390e-02,  3.33262123e-02,
       -3.09152501e-02, -5.67949461e-02, -1.36112552e-01, -2.04931791e-01,
       -1.46427808e-01,  1.65410321e-02,  2.07342753e-01,  2.58522294e-01,
        2.74758141e-01,  1.11026337e-01, -7.72118290e-02, -1.43650624e-01,
       -1.60466323e-01, -8.05993835e-02, -1.55339213e-02,  4.40076907e-02,
        3.72325816e-03,  3.23496200e-03,  4.30616169e-02,  5.52385022e-02,
        7.99890133e-02, -1.01626637e-02, -1.10538041e-01, -2.11493271e-01,
       -2.06274606e-01, -1.01626637e-01,  1.43742180e-02,  1.06234931e-01,
        6.35395367e-02, -3.36924345e-02, -1.35074923e-01, -1.45329142e-01,
       -5.12100589e-02,  7.31528672e-02,  1.78319651e-01,  1.71910764e-01,
        1.09500412e-01,  3.68053224e-02,  3.12204352e-02,  1.26499222e-01,
        2.67464217e-01,  3.86883145e-01,  4.25794244e-01,  3.60911893e-01,
       -1.09561449e-02, -2.65572069e-01, -2.87148656e-01, -3.17453536e-01,
       -6.09454634e-02,  1.44169439e-01,  1.86315500e-01,  1.38004700e-01,
        6.24103519e-02,  4.56556902e-02, -2.63069552e-02,  6.76290170e-02,
        3.98876919e-02, -1.08981597e-01, -1.59764397e-01, -2.37495041e-01,
       -1.66600543e-01, -1.38859218e-02,  1.89031648e-01,  2.54402295e-01,
        1.80791650e-01,  1.09012116e-01, -9.27152318e-02, -1.29001740e-01,
       -1.07730338e-01, -5.47196875e-02, -1.85857723e-02, -2.64290292e-02,
       -2.77413251e-02, -6.54622028e-02,  3.08542131e-02,  9.97650075e-02,
        1.25339518e-01,  8.26441237e-02, -1.24668111e-01, -2.69966735e-01,
       -3.05673391e-01, -1.85522019e-01, -6.59199805e-03,  1.34952849e-01,
        1.33152257e-01, -2.05694754e-02, -1.46641438e-01, -2.19183935e-01,
       -1.41697439e-01,  9.73540452e-03,  1.37058626e-01,  1.57383953e-01,
        9.13724174e-02,  1.88299203e-02, -2.11798456e-02,  4.38550981e-02,
        1.59825434e-01,  2.74819178e-01,  3.53373821e-01,  3.42509232e-01,
        3.00332652e-01,  1.44260994e-01, -1.31015961e-01, -2.23242897e-01,
       -2.84340953e-01, -2.36976226e-01, -2.33466597e-02,  7.13217566e-02,
        1.67363506e-01,  1.30527665e-01,  1.23203223e-01,  5.18509476e-02,
       -3.07016205e-02, -1.77922910e-02, -1.28910184e-01, -1.60161138e-01,
       -1.66386914e-01, -1.46855068e-01, -3.13119907e-02,  5.60014649e-02,
        1.78716391e-01,  1.40171514e-01,  1.12430189e-01,  2.60933256e-02,
       -8.71608631e-02, -8.75881222e-02, -1.22775964e-01, -6.97653127e-02,
       -7.46787927e-02, -1.92876980e-02, -2.04474013e-03,  2.90841395e-02,
        9.45463424e-02,  8.35596789e-02,  7.07113865e-02, -2.45674001e-02,
       -1.74016541e-01, -2.72164068e-01, -2.27332377e-01, -1.14810633e-01,
        2.41401410e-02,  1.12460707e-01,  5.88091678e-02, -7.88598285e-02,
       -1.62724693e-01, -1.93151646e-01, -9.65605640e-02,  2.38959929e-02,
        1.37089145e-01,  1.13864559e-01,  6.68355358e-02, -1.19327372e-02,
       -2.77413251e-02,  3.65001373e-02,  1.44535661e-01,  2.44056520e-01,
        2.99050874e-01,  3.05978576e-01,  2.84096805e-01,  2.17200232e-01,
       -8.36207160e-03, -1.88695944e-01, -2.44636372e-01, -3.05063021e-01,
       -1.59367656e-01,  3.52183599e-02,  1.20181890e-01,  1.86071352e-01,
        1.38584552e-01,  8.57264931e-02, -4.56556902e-02, -6.51875362e-02,
       -8.99685659e-02, -1.53660695e-01, -1.18717002e-01, -1.24423963e-01,
       -7.16879788e-02,  7.11081271e-03,  1.06265450e-01,  1.52562029e-01,
        1.14291818e-01,  8.00500504e-02, -4.14136174e-02, -1.08096561e-01,
       -1.17557299e-01, -1.08645894e-01, -5.31327250e-02, -2.58796960e-02,
        2.98165838e-02,  2.15765862e-02,  6.08233894e-02,  6.36310923e-02,
        4.38550981e-02, -4.76088748e-03, -9.28983428e-02, -2.19428083e-01,
       -2.46955779e-01, -1.62205878e-01, -5.08743553e-02,  5.61845759e-02,
        7.70287179e-02, -1.91351054e-02, -1.33030183e-01, -1.78991058e-01,
       -1.43742180e-01, -3.10983612e-02,  8.40479751e-02,  1.30588702e-01,
        7.58384960e-02,  2.01422163e-02, -2.72530290e-02,  1.15054781e-02,
        9.25016022e-02,  2.01147496e-01,  2.66884365e-01,  3.03079318e-01,
        2.97311319e-01,  2.84737693e-01,  2.01452681e-01, -6.85445723e-02,
       -2.09326456e-01, -2.90322581e-01, -2.97097690e-01, -7.63878292e-02,
        6.35395367e-02,  2.02337718e-01,  1.65746025e-01,  1.27567370e-01,
        2.60933256e-02, -8.83816034e-02, -6.21051668e-02, -1.11270486e-01,
       -1.24118778e-01, -1.25034333e-01, -1.27109592e-01, -5.21561327e-02,
        1.22074038e-02,  1.46580401e-01,  1.38889737e-01,  1.16672262e-01,
        3.79955443e-02, -8.38953825e-02, -1.10934782e-01, -1.30466628e-01,
       -5.93890194e-02, -4.04980621e-02,  1.13834040e-02,  2.14850307e-02,
        1.82500687e-02,  5.57878353e-02,  4.67238380e-02,  5.56962798e-02,
       -1.14444411e-02, -9.40885647e-02, -2.21442305e-01, -2.21991638e-01,
       -1.39469588e-01, -3.60118412e-02,  5.47196875e-02,  4.66322825e-02,
       -4.52589496e-02, -1.39591662e-01, -1.60954619e-01, -9.87273782e-02,
       -1.22074038e-04,  9.88494522e-02,  1.03030488e-01,  5.09353923e-02,
       -6.50044252e-03, -2.48420667e-02,  3.32956938e-02,  1.13742485e-01,
        2.09601123e-01,  2.46772668e-01,  2.88644063e-01,  2.71309549e-01,
        2.80526139e-01,  2.34809412e-01, -1.05594043e-02, -2.07556383e-01,
       -2.84005249e-01, -3.30332347e-01, -1.19388409e-01,  6.38142033e-02,
        2.04992828e-01,  1.95562609e-01,  1.18564409e-01,  3.94299142e-02,
       -1.01351970e-01, -7.02841273e-02, -8.15149388e-02, -1.07547227e-01,
       -1.00253304e-01, -1.20914335e-01, -5.57573168e-02,  4.54725791e-03,
        1.36539811e-01,  1.68828394e-01,  1.18350780e-01,  6.08233894e-02,
       -8.31934568e-02, -1.10171819e-01, -1.25980407e-01, -5.25528733e-02,
       -1.45268105e-02,  1.05899228e-02,  3.21665090e-02,  5.37125767e-03,
        4.37025056e-02,  4.19629505e-02,  5.71916868e-02,  2.83822138e-03,
       -7.72728660e-02, -1.77190466e-01, -2.27942747e-01, -1.59672842e-01,
       -7.56248665e-02,  2.19122898e-02,  3.45469527e-02, -2.20038453e-02,
       -1.09378338e-01, -1.47007660e-01, -1.05441450e-01, -2.45063631e-02,
        5.85345012e-02,  8.62758263e-02,  4.74868007e-02,  6.22577593e-03,
       -2.10882900e-02,  2.52998444e-02,  8.97854549e-02,  1.76976836e-01,
        2.20740379e-01,  2.61238441e-01,  2.77291177e-01,  2.86873989e-01,
        2.93282876e-01,  1.20181890e-01, -1.52897732e-01, -2.84890286e-01,
       -3.22611164e-01, -2.20221564e-01,  9.88799707e-03,  1.85827204e-01,
        2.10272530e-01,  1.31260109e-01,  5.24918363e-02, -5.35294656e-02,
       -9.32950835e-02, -5.67339091e-02, -8.00195318e-02, -1.22409742e-01,
       -1.35258034e-01, -1.02359081e-01, -2.53608814e-02,  9.21658986e-02,
        1.82073428e-01,  1.44932402e-01,  6.32038331e-02, -4.37330241e-02,
       -1.26041444e-01, -1.26560259e-01, -7.33359783e-02, -1.25736259e-02,
        8.85036775e-04,  1.26346629e-02,  7.84325694e-03,  1.70903653e-02,
        4.80056154e-02,  6.19830927e-02,  3.53404340e-02, -3.81176183e-02,
       -1.41972106e-01, -2.37464522e-01, -2.02612384e-01, -1.08890042e-01,
       -4.36414686e-03,  4.68459120e-02,  1.08951079e-02, -8.32850124e-02,
       -1.44199957e-01, -1.24057741e-01, -5.00198370e-02,  3.52488784e-02,
        8.38648640e-02,  6.04266488e-02,  1.02847377e-02, -1.62968841e-02,
        7.32444227e-03,  7.26950896e-02,  1.48075808e-01,  2.07342753e-01,
        2.27851192e-01,  2.60658589e-01,  2.72164068e-01,  2.98226875e-01,
        2.55195776e-01,  3.54014710e-03, -2.52204962e-01, -3.08206427e-01,
       -3.04818873e-01, -1.18350780e-01,  1.19052705e-01,  2.23120823e-01,
        1.81951353e-01,  7.28782006e-02,  8.63673818e-03, -8.94802698e-02,
       -7.30002747e-02, -4.12610248e-02, -1.03579821e-01, -1.35044404e-01,
       -1.34342479e-01, -5.68559832e-02,  3.84533219e-02,  1.56254769e-01,
        1.87353130e-01,  8.67946409e-02, -6.25629444e-03, -1.05929746e-01,
       -1.25156407e-01, -8.65810114e-02, -1.62663656e-02,  1.77922910e-02,
       -1.70903653e-03,  7.96533097e-03, -2.07525864e-03,  3.46385083e-02,
        5.41703543e-02,  5.56047243e-02, -7.23288675e-03, -8.13318278e-02,
       -1.77892392e-01, -2.42866298e-01, -1.56102176e-01, -5.94500565e-02,
        2.24921415e-02,  2.38959929e-02, -2.63985107e-02, -1.15665151e-01,
       -1.40232551e-01, -8.81374554e-02, -1.86468093e-02,  3.50047304e-02,
        5.66728721e-02,  2.79244362e-02, -3.54014710e-03, -3.26548051e-03,
        4.83718375e-02,  9.88494522e-02,  1.57444990e-01,  1.91137425e-01,
        2.16498306e-01,  2.56752220e-01,  2.90902432e-01,  3.16934721e-01,
        2.45307779e-01, -4.55946532e-02, -3.03964354e-01, -3.06619465e-01,
       -2.77901547e-01, -4.79750969e-02,  1.70140690e-01,  2.29071932e-01,
        1.42155217e-01,  2.10272530e-02, -7.01925718e-03, -9.04873806e-02,
       -5.41703543e-02, -3.84838404e-02, -1.29795221e-01, -1.56132694e-01,
       -1.33640553e-01, -9.88799707e-03,  8.85341960e-02,  1.87505722e-01,
        1.68401135e-01,  2.08441420e-02, -6.05792413e-02, -1.31565294e-01,
       -9.21964171e-02, -4.76699118e-02,  1.26956999e-02,  8.85036775e-03,
       -3.27158422e-02, -7.87377544e-03,  4.18103580e-03,  6.61946471e-02,
        6.46076846e-02,  4.12305063e-02, -5.53605762e-02, -1.23844111e-01,
       -1.94708090e-01, -2.00567644e-01, -1.02999969e-01, -3.05490280e-02,
        9.33866390e-03, -2.09051790e-02, -6.33869442e-02, -1.13315226e-01,
       -9.88189337e-02, -4.63576159e-02, -2.56355480e-03,  1.85552538e-02,
        2.85653249e-02,  1.33976257e-02,  1.40690329e-02,  2.72530290e-02,
        6.88192389e-02,  8.94497513e-02,  1.36234626e-01,  1.78106021e-01,
        2.37128819e-01,  2.82326731e-01,  3.25846126e-01,  3.23618274e-01,
        2.28583636e-01, -1.10415967e-01, -3.40372936e-01, -2.97952208e-01,
       -2.45033113e-01,  4.18103580e-03,  1.77251503e-01,  2.12042604e-01,
        8.75576037e-02,  9.46073794e-04,  1.65715506e-02, -6.52485733e-02,
       -3.12204352e-02, -6.81783502e-02, -1.73100986e-01, -1.96661275e-01,
       -1.20548112e-01,  4.50148015e-02,  1.24851222e-01,  1.95257424e-01,
        1.17893002e-01, -2.89925840e-02, -8.24915311e-02, -1.03610340e-01,
       -3.61033967e-02, -2.49946593e-02,  7.11081271e-03, -3.90942106e-02,
       -5.30411695e-02, -6.10370190e-04,  4.98367260e-02,  9.64995270e-02,
        4.41297647e-02, -1.85247353e-02, -1.14413892e-01, -1.38004700e-01,
       -1.71666616e-01, -1.43101291e-01, -7.04062014e-02, -3.52183599e-02,
       -2.72225105e-02, -5.30716880e-02, -7.07724235e-02, -9.68657491e-02,
       -7.85546434e-02, -4.54420606e-02, -2.06305124e-02, -1.22074038e-04,
        2.52388073e-02,  3.22580645e-02,  2.85042879e-02,  3.25327311e-02,
        5.91448714e-02,  8.48719749e-02,  1.49174474e-01,  1.94952239e-01,
        2.54554888e-01,  2.85500656e-01,  3.24289682e-01,  3.17606128e-01,
        1.83690909e-01, -2.16711936e-01, -2.95419172e-01, -2.65541551e-01,
       -1.66600543e-01,  7.39463485e-02,  1.68431654e-01,  1.70232246e-01,
       -7.72118290e-03,  4.50148015e-02,  1.56254769e-02, -3.41807306e-02,
       -2.54524369e-02, -1.24912259e-01, -2.04443495e-01, -1.96966460e-01,
       -2.27362896e-02,  9.87578967e-02,  1.40263070e-01,  1.58299509e-01,
        1.41300699e-02, -5.34073916e-02, -8.18201239e-02, -2.15460677e-02,
       -2.01422163e-02, -3.75682852e-02, -4.75783563e-02, -8.41395306e-02,
       -1.77007355e-02,  4.20850246e-02,  1.06814783e-01,  5.06302072e-02,
       -2.02032533e-02, -8.52992340e-02, -1.06692709e-01, -9.52482681e-02,
       -1.37607959e-01, -1.23355815e-01, -1.12460707e-01, -6.48823511e-02,
       -4.28479873e-02, -2.68562883e-02, -5.68254646e-02, -9.42106388e-02,
       -9.47904904e-02, -6.39362774e-02, -1.29398480e-02,  2.98165838e-02,
        4.88906522e-02,  2.92367321e-02,  6.28681295e-03,  2.82906583e-02,
        6.66219062e-02,  1.19907224e-01,  1.71330912e-01,  2.14300974e-01,
        2.42805261e-01,  2.83211768e-01,  3.21970275e-01,  3.28897977e-01,
        1.94616535e-01, -2.16467788e-01, -3.43638417e-01, -2.82753990e-01,
       -1.39774773e-01,  8.07519761e-02,  1.73528245e-01,  1.42124699e-01,
       -3.15561388e-02,  2.17291787e-02,  6.10065004e-02,  1.28788110e-02,
       -2.99081393e-02, -1.54484695e-01, -2.40394299e-01, -2.09021271e-01,
        5.37125767e-03,  1.36906034e-01,  1.37089145e-01,  1.02816858e-01,
       -1.87078463e-02, -4.97756890e-02, -3.18308054e-02,  3.41502121e-02,
       -1.19327372e-02, -8.44141972e-02, -1.07974487e-01, -8.99990844e-02,
        1.63274026e-02,  8.52076785e-02,  1.03915525e-01, -4.54725791e-03,
       -6.55842769e-02, -9.81170080e-02, -9.35087130e-02, -1.16306040e-01,
       -1.46305734e-01, -1.22135075e-01, -9.34781945e-02, -2.90536210e-02,
       -1.37943663e-02, -3.65611744e-02, -9.80559709e-02, -1.15146336e-01,
       -8.90835292e-02, -3.23496200e-02,  2.73445845e-02,  4.67238380e-02])
In [3]:
alpha = 0.95#choosing the zero for pre-emphasis
dataPreEmph = np.array([data[n] - alpha*data[n-1] for n in range(1, len(data))])#gives singal from n=1 since we start from 1
dataPreEmph = np.insert(dataPreEmph, 0, data[0])

Question 2

Compute and plot the narrowband magnitude spectrum slice using a Hamming window of duration = 30 ms on a segment near the centre of the given audio file.

In [4]:
dur = 0.03 #duration of the window in seconds
windowLength = int(sampleRate*dur) #the window length in samples (each sample is 1/sampleRate long)
totalSignalLength = int(len(dataPreEmph)) #the lenght of the entire signal in terms of samples
signalToAnalyse = dataPreEmph[totalSignalLength//2 - windowLength//2 : totalSignalLength//2 + windowLength//2] #taking the sample of the signal in the center
checkStr = "Window length matches the signal sample Proceed....." if (len(signalToAnalyse)==windowLength)  else "Check the signal sample again <Lengths do not match>"
print(checkStr) #checking whether the calculated window lenght and the lenght of windowed signal match
Window length matches the signal sample Proceed.....

Generating the hamming window

I have used the numpy implementation of the hamming window

In [5]:
window = np.hamming(windowLength)
plt.gcf().set_dpi(300)
plt.plot(window) 
plt.title("Hamming window")
plt.ylabel("Amplitude") 
plt.xlabel("Sample") 
plt.show() 

Generating the narrawband magnitude spectra

Here we multply the hamming window of appropriate length with our signal sample and take its dft. Thus, we can generate the required magnitude plot

In [6]:
dftSize = 10000 #setting the dft size
signalAfterHamming = signalToAnalyse*np.hamming(windowLength) #choosing the hamming window
signalHammingDFT = 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]) #computing the magnitudes of the signal afterwindowing and hamming
In [7]:
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
plt.plot(np.arange(0, 4000, 8000/dftSize)[0:] , signalHammingDFT)
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT"])
Out[7]:
<matplotlib.legend.Legend at 0x7fd584488110>

Question 3

With the same 30 ms segment of part 2, compute the autocorrelation coefficients required for LPC calculation at various p = 2,4,6,8,10. Use the Levinson-Durbin recursion to compute the LP coefficients from the autocorrelation coefficients. Plot error signal energy (i.e. square of gain) vs p.

For a particular p, we only need to compute the autocorrelation functions from 0 to p-1, since only these coefficients show up in our equations.

To compute the ACF :

$$ r(k) = \sum_{m = -M}^{M} x[m] \times x[m+k]$$$$ \implies r(k) = \sum_{m = 0}^{2M+1} x[m] \times x[m+k] \text{(when signal is indexed from 0)} $$

Further, since our signal is windowed i.e finite length the autocorrelation is always 0 for all $ k \geq N \text{(the window length)}$

In [8]:
def computeACF(signal, lag):
    #we assume that the signal is windowed
    #padding the signal till twice the length since the signal is finite
    signal = np.array(signal, dtype = float)
    return (np.sum(np.pad(signal, (0, 3*len(signal)), 'constant')*np.pad(signal, (lag, 3*len(signal)-lag), 'constant')))
In [9]:
k = [i for i in range(0, 240)]
x = np.linspace(0, 239, 2000)
acf = [computeACF(signalAfterHamming, i) for i in range(0, 240)]
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(k, acf,'.')
f = interp1d(k, acf, kind='cubic')
plt.plot(x, f(x), alpha = 0.3, linewidth = 0.9, color = 'g')
plt.vlines(k, 0, acf, color='r', linestyles='--', alpha = 0.3)
plt.axhline(0, color = 'k', alpha = 0.4)
plt.legend(["ACF at different lag values", "Interpolated ACF plot"])
plt.title("AutoCorrelation function at different lag values", fontsize= 15)
plt.xlabel("Lag (k)", fontsize= 15)
plt.ylabel("ACF Value", fontsize= 15)
Out[9]:
Text(0, 0.5, 'ACF Value')
In [10]:
plt.figure(figsize=(12, 16))
plt.suptitle("Signal of the sound at different positions")
plt.subplot(311)
lag = 0
plt.plot(np.pad(signalAfterHamming,(lag, 3*len(signalAfterHamming)-lag), 'constant'))
plt.xlabel("Samples[n]")
plt.ylabel("Amplitude")
plt.title("Original Signal")
plt.subplot(312)
lag = 200
plt.plot(np.pad(signalAfterHamming,(lag, 3*len(signalAfterHamming)-lag), 'constant'))
plt.xlabel("Samples[n]")
plt.ylabel("Amplitude")
plt.title(f"Shifted Signal by {lag}")
plt.subplot(313)
lag = 200
plt.plot(np.pad(signalAfterHamming,(0, 3*len(signalAfterHamming)), 'constant')*np.pad(signalAfterHamming,(lag, 3*len(signalAfterHamming)-lag), 'constant'))
plt.xlabel("Samples[n]")
plt.ylabel("Amplitude")
plt.title(f"Signal product actually obtained at {lag}")
Out[10]:
Text(0.5, 1.0, 'Signal product actually obtained at 200')

The Levinson Durban Recursion

We use the LD recursion as mentioned in class.

Emin is computed using the autocorrelation values

$$Emin =\left( r[0] - \sum_{k=1}^{P} (r[k] \times a_k) \right) \text{\;where P is the order of LD recursion used}$$$$Error = Emin = \sum_{n} \left ( s[n] - \sum_{k=1}^{P} (a_k \times s[n-k]) \right)^2$$
In [11]:
rp = np.array([computeACF(signalAfterHamming, lag = 1)])
pp = np.array([computeACF(signalAfterHamming, lag = 1)])
alpha = np.array([rp[0]/computeACF(signalAfterHamming, 0)])
beta = np.array([rp[0]/computeACF(signalAfterHamming, 0)])
coef = [alpha]

order = 40
for p in range(1, order):
    k = (computeACF(signalAfterHamming, p+1) - np.sum(pp*alpha))/(computeACF(signalAfterHamming, 0) - np.sum(rp*alpha))
    alpha = np.append(alpha, 0) - (k*np.append(beta, -1))
    beta = alpha[::-1]
    rp = np.append(rp, computeACF(signalAfterHamming, lag = p+1))
    pp = rp[::-1]
    coef.append(alpha)
In [12]:
coef
Out[12]:
[array([0.63677182]),
 array([ 0.96454776, -0.51474629]),
 array([ 0.7516907 , -0.11588805, -0.4135184 ]),
 array([ 0.68322148, -0.13507646, -0.28905557, -0.16557719]),
 array([ 0.70802401, -0.09177757, -0.26882188, -0.26791991,  0.14979436]),
 array([ 0.73630801, -0.14236591, -0.31958052, -0.28524925,  0.28348265,
        -0.18881888]),
 array([ 0.69995374, -0.08778552, -0.37450104, -0.34677974,  0.25607221,
        -0.04705368, -0.19253518]),
 array([ 0.5881573 , -0.11510746, -0.22581153, -0.54813899,  0.03861646,
        -0.09802675,  0.21389618, -0.58065459]),
 array([ 0.75507299, -0.17659432, -0.19763264, -0.55923973,  0.19618517,
        -0.03311469,  0.24698511, -0.74972702,  0.28746124]),
 array([ 0.71103761, -0.06174574, -0.2354676 , -0.55416698,  0.16613212,
         0.05255368,  0.2772599 , -0.72267503,  0.17179372,  0.15318719]),
 array([ 0.71997421, -0.05172367, -0.27762686, -0.53799225,  0.16919798,
         0.06224546,  0.24493102, -0.73641169,  0.16819161,  0.19466756,
        -0.05833779]),
 array([ 0.71734943, -0.04296503, -0.27005944, -0.5711255 ,  0.18021812,
         0.06504606,  0.25254372, -0.76061748,  0.1557004 ,  0.19234036,
        -0.02594411, -0.04499282]),
 array([ 0.71321699, -0.04534791, -0.25239364, -0.55682495,  0.11035801,
         0.08824133,  0.25851798, -0.74406506,  0.10324447,  0.16753633,
        -0.02989031,  0.02089327, -0.09184658]),
 array([ 0.71050323, -0.04473058, -0.2532768 , -0.55187481,  0.11340854,
         0.06625671,  0.26615631, -0.74145783,  0.10650518,  0.15108402,
        -0.03734769,  0.01955339, -0.07077341, -0.02954664]),
 array([ 0.71048062, -0.04478475, -0.25326183, -0.55190339,  0.11352416,
         0.06633822,  0.26558886, -0.74125413,  0.10655589,  0.15117081,
        -0.03777006,  0.01935955, -0.07080764, -0.02900288, -0.00076533]),
 array([ 0.71048463, -0.04463266, -0.25289051, -0.55200492,  0.11372223,
         0.06554547,  0.26503007, -0.73736696,  0.10516313,  0.15082293,
        -0.03836538,  0.02225376, -0.06947953, -0.02876802, -0.00449112,
         0.00524405]),
 array([ 0.71053751, -0.04467794, -0.25318056, -0.55270543,  0.1139466 ,
         0.06515866,  0.26655071, -0.73630668,  0.09772879,  0.15349504,
        -0.03770453,  0.02340034, -0.075045  , -0.03131774, -0.00494112,
         0.01240736, -0.01008228]),
 array([ 7.10541793e-01, -4.46832121e-02, -2.53178460e-01, -5.52692112e-01,
         1.13978512e-01,  6.51487120e-02,  2.66566744e-01, -7.36371944e-01,
         9.76872307e-02,  1.53808132e-01, -3.78178762e-02,  2.33726284e-02,
        -7.50934506e-02, -3.10827156e-02, -4.83346368e-03,  1.24263533e-02,
        -1.03844189e-02,  4.25218868e-04]),
 array([ 0.71062353, -0.04667942, -0.25078973, -0.55362125,  0.10800346,
         0.05071344,  0.27105968, -0.7436417 ,  0.12725388,  0.17258662,
        -0.17937121,  0.07461495, -0.06256986, -0.00917254, -0.11107788,
        -0.03624233, -0.01897391,  0.1370132 , -0.19223075]),
 array([ 6.73229943e-01, -2.00269926e-02, -2.54480623e-01, -5.60671273e-01,
         8.63960885e-02,  4.89291577e-02,  2.58888310e-01, -7.29127267e-01,
         9.23617892e-02,  2.06158947e-01, -1.54617210e-01, -7.00415753e-02,
        -9.84211638e-03,  6.92466013e-04, -9.00685646e-02, -1.43935228e-01,
        -6.77586532e-02,  1.27932911e-01, -5.39970645e-02, -1.94524498e-01]),
 array([ 0.6594864 , -0.02384199, -0.24544191, -0.56545856,  0.07622678,
         0.04256563,  0.25893723, -0.72982263,  0.08741321,  0.19523493,
        -0.14005167, -0.06351603, -0.06135641,  0.01898344, -0.08661162,
        -0.13783117, -0.10737119,  0.10995335, -0.05541201, -0.14695946,
        -0.07065199]),
 array([ 0.65033269, -0.04288215, -0.25262112, -0.55121294,  0.0623157 ,
         0.02470815,  0.24771578, -0.72736313,  0.07946384,  0.18700576,
        -0.15819685, -0.03822128, -0.0500311 , -0.0755728 , -0.05306356,
        -0.13231634, -0.09749521,  0.03669221, -0.08721161, -0.15004845,
         0.01479145, -0.12956058]),
 array([ 0.64392694, -0.04215083, -0.26003983, -0.55552486,  0.06412984,
         0.01988778,  0.24117379, -0.7299867 ,  0.07572736,  0.18453212,
        -0.16008659, -0.04604286, -0.04078515, -0.07164395, -0.08902591,
        -0.12006876, -0.09627359,  0.03977323, -0.11446472, -0.16253856,
         0.01267127, -0.09740678, -0.04944208]),
 array([ 0.64309207, -0.04379562, -0.25982586, -0.55826946,  0.06219701,
         0.02055939,  0.23954813, -0.73201416,  0.07422409,  0.18332235,
        -0.16077528, -0.04682033, -0.04348835, -0.06852796, -0.08774719,
        -0.13239521, -0.09220116,  0.04010905, -0.11338183, -0.17191907,
         0.00828027, -0.09811853, -0.03856882, -0.01688586]),
 array([ 0.64350688, -0.04284817, -0.25741554, -0.55847287,  0.06642026,
         0.02334465,  0.23856283, -0.72974921,  0.07747642,  0.18547789,
        -0.15909186, -0.04575203, -0.04233819, -0.06457846, -0.09225057,
        -0.13421855, -0.07421897,  0.03422446, -0.11388688, -0.17344696,
         0.02199437, -0.09173582, -0.03749296, -0.03268365,  0.02456536]),
 array([ 0.64603746, -0.04621505, -0.26127785, -0.56792297,  0.068686  ,
         0.00547715,  0.22683086, -0.7262236 ,  0.06983081,  0.17165147,
        -0.16859498, -0.05240453, -0.04669962, -0.06929157, -0.10863929,
        -0.11511169, -0.06623779, -0.04095009, -0.08931151, -0.17104213,
         0.0288366 , -0.14926647, -0.06401043, -0.03709762,  0.09085573,
        -0.10301423]),
 array([ 0.65092084, -0.05052205, -0.25951925, -0.56488856,  0.07576195,
         0.00411015,  0.23493909, -0.7219898 ,  0.07177205,  0.17479146,
        -0.16313813, -0.0472545 , -0.04341487, -0.06707778, -0.10615506,
        -0.10711947, -0.07437491, -0.04426041, -0.05488498, -0.18179502,
         0.02857695, -0.15252252, -0.03708811, -0.02471177,  0.09304655,
        -0.13363956,  0.04740487]),
 array([ 0.6536131 , -0.05811183, -0.25423486, -0.56629201,  0.07365561,
        -0.00455204,  0.23656206, -0.73231446,  0.06865497,  0.17227779,
        -0.1673621 , -0.05333812, -0.04944372, -0.07088732, -0.10862072,
        -0.10980319, -0.08363999, -0.0343335 , -0.05080884, -0.2227989 ,
         0.04191982, -0.15228909, -0.03278537, -0.05679342,  0.0783077 ,
        -0.13650885,  0.08437254, -0.05679288]),
 array([ 0.65278817, -0.05688629, -0.25621769, -0.56515457,  0.07283067,
        -0.00502826,  0.23435002, -0.73170557,  0.06541876,  0.17153978,
        -0.1678608 , -0.05455301, -0.05103864, -0.07246506, -0.10965037,
        -0.11052137, -0.08441474, -0.03676448, -0.04830646, -0.22180167,
         0.03128278, -0.14885297, -0.03285149, -0.05572356,  0.07008217,
        -0.14020167,  0.08352845, -0.047299  , -0.01452524]),
 array([ 0.65204805, -0.05929635, -0.25196162, -0.57229835,  0.07640161,
        -0.00786757,  0.23267612, -0.73929016,  0.06701273,  0.16023819,
        -0.17032218, -0.05642629, -0.05533987, -0.07809652, -0.11523745,
        -0.11421372, -0.08701534, -0.03954415, -0.05685956, -0.21306111,
         0.0346161 , -0.18613597, -0.02091053, -0.05597977,  0.07379316,
        -0.16899831,  0.07047325, -0.05019755,  0.01873664, -0.05095356]),
 array([ 0.65316132, -0.05970572, -0.25086487, -0.57383809,  0.080094  ,
        -0.00947985,  0.23389921, -0.73883329,  0.07107955,  0.15948188,
        -0.16566709, -0.05518399, -0.05447588, -0.07619535, -0.11274204,
        -0.11169594, -0.08530904, -0.03833504, -0.05562673, -0.2093398 ,
         0.03111511, -0.18760011, -0.00475804, -0.06106343,  0.07396505,
        -0.17066758,  0.0829772 , -0.04469253,  0.02003218, -0.06519993,
         0.02184865]),
 array([ 6.52416688e-01, -5.74836150e-02, -2.51547591e-01, -5.72314911e-01,
         7.72660195e-02, -3.66326904e-03,  2.31378377e-01, -7.36752165e-01,
         7.12417106e-02,  1.65875544e-01, -1.66727531e-01, -4.80494018e-02,
        -5.25800478e-02, -7.48888415e-02, -1.09834588e-01, -1.07889189e-01,
        -8.14666375e-02, -3.57382027e-02, -5.37701145e-02, -2.07459056e-01,
         3.67612679e-02, -1.93035470e-01, -7.18052309e-03, -3.58829843e-02,
         6.59934492e-02, -1.70344493e-01,  8.02474863e-02, -2.51353528e-02,
         2.85819969e-02, -6.31650737e-02, -4.11974389e-04,  3.40813569e-02]),
 array([ 6.52597688e-01, -5.74858029e-02, -2.51883050e-01, -5.72163117e-01,
         7.71325301e-02, -3.23708870e-03,  2.30473707e-01, -7.36401686e-01,
         7.10511424e-02,  1.65837409e-01, -1.67752709e-01, -4.78541691e-02,
        -5.36818265e-02, -7.51744052e-02, -1.10024387e-01, -1.08321844e-01,
        -8.20396181e-02, -3.63215150e-02, -5.41678360e-02, -2.07738300e-01,
         3.65060860e-02, -1.93920931e-01, -6.29958714e-03, -3.55046320e-02,
         6.20806875e-02, -1.69115683e-01,  8.02280313e-02, -2.47250065e-02,
         2.55425327e-02, -6.45009989e-02, -7.17259795e-04,  3.75462277e-02,
        -5.31082477e-03]),
 array([ 0.65177114, -0.05164233, -0.25199468, -0.58220168,  0.08110782,
        -0.00708515,  0.24295993, -0.76272187,  0.08071302,  0.16031167,
        -0.16873314, -0.0780349 , -0.04800023, -0.10750559, -0.11845476,
        -0.11397471, -0.09480779, -0.0531801 , -0.07129139, -0.21943801,
         0.02815136, -0.20136868, -0.03240765, -0.00969466,  0.07313868,
        -0.28372498,  0.11609762, -0.02522881,  0.03754699, -0.15354915,
        -0.03991888,  0.02859947,  0.0962557 , -0.15563421]),
 array([ 0.63041252, -0.03843257, -0.2480698 , -0.58767999,  0.06003534,
        -0.00193235,  0.23949763, -0.74678909,  0.04177573,  0.17034894,
        -0.1700636 , -0.0824824 , -0.07563527, -0.10364221, -0.14856956,
        -0.12375846, -0.10210602, -0.06619114, -0.08693283, -0.23569427,
         0.01339772, -0.20795604, -0.04311685, -0.03285093,  0.09513921,
        -0.27264824,  0.0114247 ,  0.00811405,  0.03657466, -0.14241824,
        -0.11981793, -0.00598328,  0.08916851, -0.06618772, -0.13723604]),
 array([ 0.62955763, -0.03884488, -0.24751434, -0.58771726,  0.05928896,
        -0.00281952,  0.23972547, -0.74673855,  0.0418469 ,  0.16865052,
        -0.16947094, -0.08268704, -0.07590386, -0.10493764, -0.1484861 ,
        -0.12522668, -0.10264755, -0.06660347, -0.08756889, -0.2364652 ,
         0.01247223, -0.20860166, -0.04358801, -0.03336474,  0.09407983,
        -0.27158708,  0.01168493,  0.00346205,  0.03806657, -0.14243028,
        -0.11944395, -0.00964414,  0.0876232 , -0.06642713, -0.13330898,
        -0.00622934]),
 array([ 0.62976822, -0.03433827, -0.24526872, -0.59067943,  0.05961499,
         0.00121838,  0.24454043, -0.74802542,  0.04172986,  0.1682555 ,
        -0.16028973, -0.08586748, -0.07477594, -0.10346411, -0.14143416,
        -0.12564832, -0.09465366, -0.06364314, -0.0853173 , -0.23299512,
         0.01670561, -0.20358197, -0.04004051, -0.03079875,  0.09687513,
        -0.26585799,  0.00598357,  0.00204738,  0.06331064, -0.15053438,
        -0.11934863, -0.01164845,  0.10749144, -0.05805972, -0.1319958 ,
        -0.02751202,  0.03380577]),
 array([ 0.62380655, -0.02948651, -0.22199118, -0.58044056,  0.04065881,
         0.00327259,  0.26558763, -0.72147859,  0.03056499,  0.16789444,
        -0.16134494, -0.03898327, -0.09185992, -0.09803274, -0.13437299,
        -0.08974652, -0.09759971, -0.02255431, -0.07027155, -0.22177162,
         0.03339784, -0.18142382, -0.01509851, -0.01255279,  0.11006191,
        -0.25071521,  0.03425076, -0.02762457,  0.05595155, -0.01861966,
        -0.16247348, -0.01186331,  0.0969783 ,  0.04610694, -0.08874252,
        -0.02145644, -0.07725423,  0.17635058]),
 array([ 0.62898191, -0.03175369, -0.22262087, -0.58304489,  0.04201191,
         0.00611861,  0.26523948, -0.7262467 ,  0.03001856,  0.16953645,
        -0.16215564, -0.03797811, -0.09921766, -0.09480275, -0.13474138,
        -0.09018961, -0.10292395, -0.02157419, -0.07677987, -0.22383388,
         0.03273594, -0.18428808, -0.0177323 , -0.01649623,  0.10718495,
        -0.25341102,  0.03310672, -0.03235956,  0.06087875, -0.01772267,
        -0.1836467 , -0.00406911,  0.09707434,  0.04730016, -0.1057767 ,
        -0.02797121, -0.07811956,  0.19465742, -0.02934698]),
 array([ 0.6316245 , -0.0492819 , -0.21558648, -0.58052618,  0.05153672,
         0.0018594 ,  0.25649828, -0.72588029,  0.0465553 ,  0.17113232,
        -0.16763755, -0.03506425, -0.1021988 , -0.07198399, -0.144393  ,
        -0.08870418, -0.10132722, -0.0049797 , -0.07972763, -0.20367843,
         0.03964969, -0.1823454 , -0.00846436, -0.00837498,  0.11931793,
        -0.24487437,  0.04204092, -0.02893976,  0.07548029, -0.03298883,
        -0.18634976,  0.06132683,  0.07319046,  0.0467492 , -0.10955972,
         0.02452991, -0.05807334,  0.19751673, -0.08598457,  0.09004646])]
In [13]:
#checking whether the coefficient are correct
#should preserve the acf!
for p in range(1, len(coef)+1):
    computeACF(signalAfterHamming, 1) == coef[0][0]*computeACF(signalAfterHamming, 0)
    Rind = np.array([[ np.abs(i-k) for k in range(1, p+1)] for i in range(1, p+1)])
    R = np.array([[computeACF(signalAfterHamming, np.abs(i-k)) for k in range(1, p+1)] for i in range(1, p+1)])
    acfArr = np.array([computeACF(signalAfterHamming, i) for i in range(1, p+1)])
    coefArr = coef[p-1]
    if np.sum((np.round(np.matmul(R, coefArr) - acfArr, 7))**2)==0:
        print(f"Order {p} coefficients verified and correct...proceed")
    else:
        print("Error <Coefficients don't match>")
        break
Order 1 coefficients verified and correct...proceed
Order 2 coefficients verified and correct...proceed
Order 3 coefficients verified and correct...proceed
Order 4 coefficients verified and correct...proceed
Order 5 coefficients verified and correct...proceed
Order 6 coefficients verified and correct...proceed
Order 7 coefficients verified and correct...proceed
Order 8 coefficients verified and correct...proceed
Order 9 coefficients verified and correct...proceed
Order 10 coefficients verified and correct...proceed
Order 11 coefficients verified and correct...proceed
Order 12 coefficients verified and correct...proceed
Order 13 coefficients verified and correct...proceed
Order 14 coefficients verified and correct...proceed
Order 15 coefficients verified and correct...proceed
Order 16 coefficients verified and correct...proceed
Order 17 coefficients verified and correct...proceed
Order 18 coefficients verified and correct...proceed
Order 19 coefficients verified and correct...proceed
Order 20 coefficients verified and correct...proceed
Order 21 coefficients verified and correct...proceed
Order 22 coefficients verified and correct...proceed
Order 23 coefficients verified and correct...proceed
Order 24 coefficients verified and correct...proceed
Order 25 coefficients verified and correct...proceed
Order 26 coefficients verified and correct...proceed
Order 27 coefficients verified and correct...proceed
Order 28 coefficients verified and correct...proceed
Order 29 coefficients verified and correct...proceed
Order 30 coefficients verified and correct...proceed
Order 31 coefficients verified and correct...proceed
Order 32 coefficients verified and correct...proceed
Order 33 coefficients verified and correct...proceed
Order 34 coefficients verified and correct...proceed
Order 35 coefficients verified and correct...proceed
Order 36 coefficients verified and correct...proceed
Order 37 coefficients verified and correct...proceed
Order 38 coefficients verified and correct...proceed
Order 39 coefficients verified and correct...proceed
Order 40 coefficients verified and correct...proceed
In [14]:
#computing the minima energy directly from formula
Emin = []
for p in range(1,len(coef)+1):
    acfArr = np.array([computeACF(signalAfterHamming, i) for i in range(1, p+1)])
    E = computeACF(signalAfterHamming, 0) - np.sum(acfArr*coef[p-1])
    Emin.append(E)
In [15]:
error = []
for p in range(1, len(coef)+1):
    signalShifted = np.pad(np.array(signalAfterHamming, dtype = float),(p, 3*len(signalAfterHamming)-p), 'constant')
    signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(0, 3*len(signalAfterHamming)), 'constant')
    signalEstimate = np.zeros_like(signalShifted)
    for i in range(0, len(signalAfterHamming) + p-1):
        signalEstimate[i] = np.sum(signalShifted[i:(i+p)]*coef[p-1][::-1])
    error.append(np.sum((signalToEstimate - signalEstimate)**2))

Error plots

In [16]:
plt.figure(figsize=(16, 8))
plt.gcf().set_dpi(300)
plt.suptitle("Error computed for different orders of LPA for different methods of error computation")
plt.subplot(121)
plt.plot([i for i in range(1, order + 1)], error, 'ro')
plt.plot([i for i in range(1, order + 1)],  error, 'b', alpha = 0.3, linewidth = 1)
plt.vlines([i for i in range(1, order + 1)] , 0, error, color='g', linestyles='--', alpha = 0.3)
plt.ylim(bottom=0)
plt.title("Calculated MSE error/loss from the computed coefficients")
plt.xlabel("Order of LPA (p)")
plt.ylabel("Calculated error")
plt.legend(["Error computed"])
plt.subplot(122)
plt.plot([i for i in range(1, order + 1)], Emin, 'ro')
plt.plot([i for i in range(1, order + 1)], Emin, 'b', alpha = 0.3, linewidth = 1)
plt.vlines([i for i in range(1, order + 1)] , 0, Emin, color='g', linestyles='--', alpha = 0.3)
plt.ylim(bottom=0)
plt.title("Calculated MSE error/loss from direct loss equation")
plt.xlabel("Order of LPA (p)")
plt.ylabel("Cmputed error from formula")
plt.legend(["Error computed from formula"])
Out[16]:
<matplotlib.legend.Legend at 0x7fd581988190>

Inference

  1. We see that the computed loss reduce as p increases and finally settle at minima around p=40
  2. This is expected since higher order lpa approximated the system better than lower order lpa

Question 4

Show the pole-zero plots of the estimated all-pole filter for p=6,10; comment.

Order of lpa p=6

In [17]:
z, p, k =tf2zpk(np.sqrt(Emin[5]), np.insert((-1)*coef[5] , 0, 1))
In [18]:
#Zeros
z
Out[18]:
array([], dtype=float64)
In [19]:
#poles
p
Out[19]:
array([-0.63742054+0.47891212j, -0.63742054-0.47891212j,
        0.65982201+0.58240455j,  0.65982201-0.58240455j,
        0.34575253+0.51376466j,  0.34575253-0.51376466j])
In [20]:
angle = np.linspace( 0 , 2 * np.pi , 150 ) 

radius = 1

x = radius * np.cos( angle ) 
y = radius * np.sin( angle )
plt.figure(figsize=(10,10)) 
plt.plot(x, y, 'r', alpha = 0.2)
plt.plot([i.real for i in p], [j.imag for j in p], 'x', label= "poles")
plt.xlim([-1.1, 1.1])
plt.ylim([-1.1, 1.1])
plt.title("Pole zero plot for LPA of order 6")
plt.xlabel("Real axis")
plt.ylabel("Imaginary axis")
plt.legend()
Out[20]:
<matplotlib.legend.Legend at 0x7fd5818cfd50>

Order of lpa p=10

In [21]:
z, p, k =tf2zpk(np.sqrt(Emin[9]),np.insert((-1)*coef[9] , 0, 1))
In [22]:
#Zeros
z
Out[22]:
array([], dtype=float64)
In [23]:
#poles
p
Out[23]:
array([ 0.56120875+0.78976252j,  0.56120875-0.78976252j,
        0.80425681+0.52217128j,  0.80425681-0.52217128j,
        0.66327132+0.j        , -0.33552838+0.87509787j,
       -0.33552838-0.87509787j, -0.83638979+0.44522166j,
       -0.83638979-0.44522166j, -0.33932847+0.j        ])
In [24]:
angle = np.linspace( 0 , 2 * np.pi , 150 ) 

radius = 1

x = radius * np.cos( angle ) 
y = radius * np.sin( angle )
plt.figure(figsize=(10,10)) 
plt.plot(x, y, 'r', alpha = 0.2)
plt.plot([i.real for i in p], [j.imag for j in p], 'x', label = "Poles")
plt.xlim([-1.1, 1.1])
plt.ylim([-1.1, 1.1])
plt.title("Pole zero plot for LPA of order 10")
plt.xlabel("Real axis")
plt.ylabel("Imaginary axis")
plt.legend()
Out[24]:
<matplotlib.legend.Legend at 0x7fd5843a1210>

Comments:

  1. We see that most poles lie very close to the unit circle in both the cases(for p = 6 a few poles are a little inside the unit circle)
  2. We see 2 poles lying on the real axis in the case of p=10 which is unexpected!

Question 5

  1. Compute the gain
    and
  2. plot the LPC spectrum magnitude (i.e. the dB magnitude frequency response of the estimated all-pole filter) for each order "p".
  3. Comment on the characteristics of the spectral envelope estimates.
  4. Comment on their shapes with reference to the short-time magnitude spectrum computed in part 2.

Gain

Gain is obtained as follows $G = \sqrt{E_{min}}$

In [25]:
#can we obtained by sqrt(Emin)

for i in range(0, order):
    print(f"Gain for order {i+1} lpa is {np.sqrt(Emin[i])}")
Gain for order 1 lpa is 0.6121560110899257
Gain for order 2 lpa is 0.5248273975246246
Gain for order 3 lpa is 0.4778531209543422
Gain for order 4 lpa is 0.471257234668652
Gain for order 5 lpa is 0.46594012146053576
Gain for order 6 lpa is 0.45755875829613224
Gain for order 7 lpa is 0.4489978668853352
Gain for order 8 lpa is 0.36555161851205675
Gain for order 9 lpa is 0.3501225154561596
Gain for order 10 lpa is 0.345990085411915
Gain for order 11 lpa is 0.34540083003193434
Gain for order 12 lpa is 0.3450510460965086
Gain for order 13 lpa is 0.3435925739538078
Gain for order 14 lpa is 0.3434425623362686
Gain for order 15 lpa is 0.34344246175445886
Gain for order 16 lpa is 0.34343773938352806
Gain for order 17 lpa is 0.3434202832925263
Gain for order 18 lpa is 0.3434202522454279
Gain for order 19 lpa is 0.33701537986354735
Gain for order 20 lpa is 0.3305775975491967
Gain for order 21 lpa is 0.3297514926536113
Gain for order 22 lpa is 0.32697218496578684
Gain for order 23 lpa is 0.326572295564563
Gain for order 24 lpa is 0.3265257341187591
Gain for order 25 lpa is 0.3264271971225869
Gain for order 26 lpa is 0.32469056593986734
Gain for order 27 lpa is 0.32432553488576416
Gain for order 28 lpa is 0.32380206751669344
Gain for order 29 lpa is 0.3237679074059267
Gain for order 30 lpa is 0.32334734062955356
Gain for order 31 lpa is 0.32327015430887424
Gain for order 32 lpa is 0.3230823543304585
Gain for order 33 lpa is 0.3230777980520862
Gain for order 34 lpa is 0.3191410163072854
Gain for order 35 lpa is 0.3161214233231459
Gain for order 36 lpa is 0.3161152897736016
Gain for order 37 lpa is 0.31593460513566596
Gain for order 38 lpa is 0.31098309512729017
Gain for order 39 lpa is 0.310849149963949
Gain for order 40 lpa is 0.30958634584683253

LPC spectra magnitude plots

For order p = 2

In [26]:
p = 2

signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
    signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])

w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]),  0, 1), None, 1)

plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))  

plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[26]:
<matplotlib.legend.Legend at 0x7fd58544e3d0>

For order p = 4

In [27]:
p = 4

signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
    signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])

w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]),  0, 1), None, 1)

plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))  

plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[27]:
<matplotlib.legend.Legend at 0x7fd581785dd0>

For order p = 6

In [28]:
p = 6

signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
    signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])

w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]),  0, 1), None, 1)

plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))  

plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[28]:
<matplotlib.legend.Legend at 0x7fd58171cfd0>

For order p = 8

In [29]:
p = 8

signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
    signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])

w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]),  0, 1), None, 1)

plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))  

plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[29]:
<matplotlib.legend.Legend at 0x7fd58169de50>

For order p = 10

In [30]:
p = 10

signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
    signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])

w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]),  0, 1), None, 1)

plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))  

plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[30]:
<matplotlib.legend.Legend at 0x7fd5815bb690>

We clearly see that as the order increases of the lpa analysis we get our estimated envelope more and more closer to the actual waveform spectra

This means that we are able to better identfy the formants (since hidden envelope is estimated) as the order of the lpa increases

Question 6

Based on the 10th-order LP coefficients, carry out the inverse filtering of the /a/ vowel segment to obtain the residual error signal. Can you measure the pitch period of the voiced sound from the residual waveform? Use the acf to detect the pitch. Compare the acf plots of the original speech and residual signals.

In [31]:
p = 10
signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(p, 3*len(signalAfterHamming)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalAfterHamming) + p-1):
    signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])
In [32]:
plt.figure(figsize = (16, 8))
plt.gcf().set_dpi(300)
plt.suptitle(f"Plots of the signal, its estimate and the error obtained in time domain for order {p}")
plt.subplot(121)
plt.plot(signalToEstimate[p:240+p], label = "Signal(Ground truth)")
plt.plot(signalEstimate[0:240], label = "Signal estimated")
plt.title("Signal along with its estimate")
plt.xlabel("Time(sec)")
plt.ylabel("Amplitude")
plt.legend()

plt.subplot(122)
plt.plot((signalEstimate -signalToEstimate)[0:240], label = "Error") 
plt.title("Error in the estimate")
plt.xlabel("Time(sec)")
plt.ylabel("Amplitude")
plt.legend()
plt.legend()
Out[32]:
<matplotlib.legend.Legend at 0x7fd5815a7790>
In [33]:
k = [i for i in range(0, 300)]
x = np.linspace(0, 299, 2000)
acf = [computeACF(signalEstimate -signalToEstimate, i) for i in range(0, 300)]
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(k, acf,'.')
f = interp1d(k, acf, kind='cubic')
plt.plot(x, f(x), alpha = 0.3, linewidth = 0.9, color = 'g')
plt.vlines(k, 0, acf, color='r', linestyles='--', alpha = 0.3)
plt.axhline(0, color = 'k', alpha = 0.4)
plt.legend(["ACF at different lag values", "Interpolated ACF plot"])
plt.title("AutoCorrelation function at different lag values of the error signal", fontsize= 15)
plt.xlabel("Lag (k)", fontsize= 15)
plt.ylabel("ACF Value", fontsize= 15)
Out[33]:
Text(0, 0.5, 'ACF Value')

Referring the Question 3 we have the original ACF plot of the signal.

Comparing the two, we have :

  1. Both of them peak at nearly 59 samples (i.e 59 $\times \frac{1}{8000} = 7.375$ micro-seconds).
  2. The above behaviour is expected since, both signals have the underlying pitch(glottal vabrations)
  3. Thus, the estimated pitch is $\frac{1}{7.375ms} = 135.925Hz$
  4. The nature of the two ACF graphs is different (i.e.) their shape etc and this is obvious since the underlying signal is different but still point (1) holds

Question 7

LP re-synthesis: We analysed a natural speech sound /a/ above.

Using a suitable set of parameter estimates as obtained there, we wish to reconstruct the sound. That is, use the best estimated LP filter with an ideal impulse train of the estimated pitch period as source excitation.

Carry out de-emphasis on the output waveform. Set the duration of the synthesized sound to be 300 ms at 8 kHz sampling frequency and view the waveform as well as listen to your created sound.

Comment on the similarity with the original sound. Try out voice modification using this analysis-synthesis method (e.g. change the voice pitch).

We will use the order 10 filter.

Say u[n] is our impulse train. Clearly, u[n] = 1 when n = $k \times T$ else u[n] = 0 where T is the pitch of the signal i.e. 59 samples as calculated in the previous question.

To calculate the output(estimated signal), we need to set up the difference equation as follows:

We know our filter looks like : $$A(z) = \frac{G}{1 - \sum_{k = 1}^{k=P}a_k z^{-k}}$$ $$\implies y[n] - \sum_{k = 1}^{k=P}a_k y[n-k] = Gx[n]$$

$$\implies y[n] = Gx[n] + \sum_{k = 1}^{k=P}a_k y[n-k]$$

We assume that the signal y[n] is causal and x[n] is our impulse i.e. impulse train

In [34]:
signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(0, 10*len(signalAfterHamming)), 'constant')
signalReconstructed = np.zeros_like(signalToEstimate)
inputImpulseTrain = np.zeros_like(signalToEstimate)
inputImpulseTrain[[i for i in range(0, signalToEstimate.shape[0], 59)]] = 1
#inputImpulseTrain[0]=1
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(inputImpulseTrain, 'o')
plt.vlines(np.arange(0, signalToEstimate.shape[0], 59), 0, 1, colors = 'r', linestyles='--')
plt.xlabel("Times(samples)")
plt.ylabel("Amplitude")
plt.title("Impulse train generated")
Out[34]:
Text(0.5, 1.0, 'Impulse train generated')
In [35]:
P = 10#choosing the order of analysis
G = np.sqrt(Emin[P-1])
inputImpulseTrain = np.pad(np.array(inputImpulseTrain, dtype = float),(P, P), 'constant', constant_values=0)
#signal y[n] is assumed casual
for i in range(P, signalToEstimate.shape[0]):
    signalReconstructed[i] = G*inputImpulseTrain[i] + np.sum(signalReconstructed[i-P:i]*coef[P-1][::-1])

We now need to de-emphasis, we can do this by passing the signal through a filter : $$H(z) = \frac{1}{(1 - \alpha z^{-1})} \text{ with $\alpha$ very close to 1, typically = 0.95}$$ $$\implies s_{de-emph}[n] = s_{emph}[n] + \alpha s_{de-emph}[n-1]$$

In [36]:
alpha = 0.95#setting the same value as pre-emphasis
signalReconstructed = signalReconstructed[P:]
signalReconstructedDeEmph = np.pad(np.zeros_like(signalReconstructed), (1, 0), 'constant')
for j in range(1, signalReconstructedDeEmph.shape[0]-1):
    signalReconstructedDeEmph[j] = signalReconstructed[j] + alpha*signalReconstructedDeEmph[j-1]
In [37]:
plt.figure(figsize=(16,8))
plt.gcf().set_dpi(300)
plt.suptitle("Signal estimated along with ground truth")
plt.subplot(121)
plt.plot([i*(1.0/8000) for i in range(len(signalReconstructed))][0:240],signalReconstructedDeEmph[1:][0:240])
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Estimate of the signal(after resynthesis and de-emphasis)")
plt.subplot(122)
plt.plot([i*(1.0/8000) for i in range(len(signalReconstructed))][0:240],signalToAnalyse)
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Ground truth(true signal")
Out[37]:
Text(0.5, 1.0, 'Ground truth(true signal')
In [38]:
scaled = np.int16(signalReconstructedDeEmph/np.max(np.abs(signalReconstructedDeEmph)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('aaReconstructed.wav', 8000, scaled)

We see the sound is very similar to the original sound

Further the waveforms are similar looking -notice their amplitudes etc. may not match exaclty

Changing pitch of the sound

Notice that the beauty of the lpa analysis is that we are finding the envelope of the sound that is the glottal vibrations are not taken into account. This means that all information about the glottal vibrations that is the pitch is encoded in the e[n] signal that is fed into our model. Here e[n] is the ideal impulse train where the distance(in samples) between two impulses encodes the information about the pitch. Thus, modifying the impulse train basically gives us a rough estimate of the sound at a different pitch

Changing pitch to higher pitch

We would need to reduce the T0 between each impulse thus F0 increases and so does Pitch

In [39]:
signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(0, 10*len(signalAfterHamming)), 'constant')
signalReconstructed = np.zeros_like(signalToEstimate)
inputImpulseTrain = np.zeros_like(signalToEstimate)
inputImpulseTrain[[i for i in range(0, signalToEstimate.shape[0], 20)]] = 1
#inputImpulseTrain[0]=1
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(inputImpulseTrain, 'o')
plt.vlines(np.arange(0, signalToEstimate.shape[0], 20), 0, 1, colors = 'r', linestyles='--')
plt.xlabel("Times(samples)")
plt.ylabel("Amplitude")
plt.title("Impulse train generated")
Out[39]:
Text(0.5, 1.0, 'Impulse train generated')
In [40]:
P = 10#choosing the order of analysis
G = np.sqrt(Emin[P-1])
inputImpulseTrain = np.pad(np.array(inputImpulseTrain, dtype = float),(P, P), 'constant', constant_values=0)
#signal y[n] is assumed casual
for i in range(P, signalToEstimate.shape[0]):
    signalReconstructed[i] = G*inputImpulseTrain[i] + np.sum(signalReconstructed[i-P:i]*coef[P-1][::-1])
In [41]:
alpha = 0.95#setting the same value as pre-emphasis
signalReconstructed = signalReconstructed[P:]
signalReconstructedDeEmph = np.pad(np.zeros_like(signalReconstructed), (1, 0), 'constant')
for j in range(1, signalReconstructedDeEmph.shape[0]-1):
    signalReconstructedDeEmph[j] = signalReconstructed[j] + alpha*signalReconstructedDeEmph[j-1]
In [42]:
plt.figure(figsize=(16,8))
plt.gcf().set_dpi(300)
plt.suptitle("Signal estimated along with ground truth")
plt.subplot(121)
plt.plot([i*(1.0/8000) for i in range(len(signalReconstructed))][0:240],signalReconstructedDeEmph[1:][0:240])
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Estimate of the signal(after resynthesis and de-emphasis)")
plt.subplot(122)
plt.plot([i*(1.0/8000) for i in range(len(signalReconstructed))][0:240],signalToAnalyse)
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Estimate of the signal(after resynthesis)")
Out[42]:
Text(0.5, 1.0, 'Estimate of the signal(after resynthesis)')
In [43]:
scaled = np.int16(signalReconstructedDeEmph/np.max(np.abs(signalReconstructedDeEmph)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('aaReconstructedPitchVariationHigh.wav', 8000, scaled)

Changing pitch to lower pitch

We would need to increase the T0 between each impulse thus F0 decreases and so does Pitch

In [44]:
signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(0, 10*len(signalAfterHamming)), 'constant')
signalReconstructed = np.zeros_like(signalToEstimate)
inputImpulseTrain = np.zeros_like(signalToEstimate)
inputImpulseTrain[[i for i in range(0, signalToEstimate.shape[0], 100)]] = 1
#inputImpulseTrain[0]=1
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(inputImpulseTrain, 'o')
plt.vlines(np.arange(0, signalToEstimate.shape[0], 100), 0, 1, colors = 'r', linestyles='--')
plt.xlabel("Times(samples)")
plt.ylabel("Amplitude")
plt.title("Impulse train generated")
Out[44]:
Text(0.5, 1.0, 'Impulse train generated')
In [45]:
P = 10#choosing the order of analysis
G = np.sqrt(Emin[P-1])
inputImpulseTrain = np.pad(np.array(inputImpulseTrain, dtype = float),(P, P), 'constant', constant_values=0)
#signal y[n] is assumed casual
for i in range(P, signalToEstimate.shape[0]):
    signalReconstructed[i] = G*inputImpulseTrain[i] + np.sum(signalReconstructed[i-P:i]*coef[P-1][::-1])
In [46]:
alpha = 0.95#setting the same value as pre-emphasis
signalReconstructed = signalReconstructed[P:]
signalReconstructedDeEmph = np.pad(np.zeros_like(signalReconstructed), (1, 0), 'constant')
for j in range(1, signalReconstructedDeEmph.shape[0]-1):
    signalReconstructedDeEmph[j] = signalReconstructed[j] + alpha*signalReconstructedDeEmph[j-1]
In [47]:
plt.figure(figsize=(16,8))
plt.gcf().set_dpi(300)
plt.suptitle("Signal estimated along with ground truth(shown for few seconds for better clarity)")
plt.subplot(121)
plt.plot([i*(1.0/8000) for i in range(len(signalReconstructed))][0:240],signalReconstructedDeEmph[1:][0:240])
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Estimate of the signal(after resynthesis and de-emphasis)")
plt.subplot(122)
plt.plot([i*(1.0/8000) for i in range(len(signalReconstructed))][0:240],signalToAnalyse)
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Estimate of the signal(after resynthesis)")
Out[47]:
Text(0.5, 1.0, 'Estimate of the signal(after resynthesis)')
In [48]:
scaled = np.int16(signalReconstructedDeEmph/np.max(np.abs(signalReconstructedDeEmph)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('aaReconstructedPitchVariationLow.wav', 8000, scaled)

[Bonus]

Perform LP analysis on the provided /s/ sound (sampled at 16 kHz).

Same steps are repeated as above

In [49]:
pathToaaSound = r"ss.wav" #path to the "pani" sound of machali
sampleRate, data = wavfile.read(pathToaaSound) #storing the used sample rate and the wav file in variables for analysis
data = data/(32767.0)
data
Out[49]:
array([-0.0212714 , -0.03814814, -0.01904355, ...,  0.04434339,
        0.03405866,  0.02053896])
In [50]:
dur = 0.03 #duration of the window in seconds
windowLength = int(sampleRate*dur) #the window length in samples (each sample is 1/sampleRate long)
totalSignalLength = int(len(data)) #the lenght of the entire signal in terms of samples
signalToAnalyse = data[totalSignalLength//2 - windowLength//2 : totalSignalLength//2 + windowLength//2] #taking the sample of the signal in the center
checkStr = "Window length matches the signal sample Proceed....." if (len(signalToAnalyse)==windowLength)  else "Check the signal sample again <Lengths do not match>"
print(checkStr) #checking whether the calculated window lenght and the lenght of windowed signal match
Window length matches the signal sample Proceed.....
In [51]:
dftSize = 10000 #setting the dft size
signalAfterHamming = signalToAnalyse*np.hamming(windowLength) #choosing the hamming window
signalHammingDFT = 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]) #computing the magnitudes of the signal afterwindowing and hamming
In [52]:
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
plt.plot(np.arange(0, 4000, 8000/dftSize)[0:] , signalHammingDFT)
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT"])
Out[52]:
<matplotlib.legend.Legend at 0x7fd581889c50>
In [53]:
k = [i for i in range(0, 240)]
x = np.linspace(0, 239, 2000)
acf = [computeACF(signalAfterHamming, i) for i in range(0, 240)]
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(k, acf,'.')
f = interp1d(k, acf, kind='cubic')
plt.plot(x, f(x), alpha = 0.3, linewidth = 0.9, color = 'g')
plt.vlines(k, 0, acf, color='r', linestyles='--', alpha = 0.3)
plt.axhline(0, color = 'k', alpha = 0.4)
plt.legend(["ACF at different lag values", "Interpolated ACF plot"])
plt.title("AutoCorrelation function at different lag values", fontsize= 15)
plt.xlabel("Lag (k)", fontsize= 15)
plt.ylabel("ACF Value", fontsize= 15)
Out[53]:
Text(0, 0.5, 'ACF Value')
In [54]:
rp = np.array([computeACF(signalAfterHamming, lag = 1)])
pp = np.array([computeACF(signalAfterHamming, lag = 1)])
alpha = np.array([rp[0]/computeACF(signalAfterHamming, 0)])
beta = np.array([rp[0]/computeACF(signalAfterHamming, 0)])
coef = [alpha]

order = 40
for p in range(1, order):
    k = (computeACF(signalAfterHamming, p+1) - np.sum(pp*alpha))/(computeACF(signalAfterHamming, 0) - np.sum(rp*alpha))
    alpha = np.append(alpha, 0) - (k*np.append(beta, -1))
    beta = alpha[::-1]
    rp = np.append(rp, computeACF(signalAfterHamming, lag = p+1))
    pp = rp[::-1]
    coef.append(alpha)
In [55]:
coef
Out[55]:
[array([0.02798905]),
 array([ 0.04839422, -0.72904164]),
 array([ 0.23663049, -0.74153688,  0.25819687]),
 array([ 0.21178449, -0.67017959,  0.23542618,  0.09622891]),
 array([ 0.13741194, -0.8521337 ,  0.75338867, -0.06745322,  0.77287119]),
 array([-0.03786684, -0.83683604,  0.58252832,  0.12580143,  0.74170766,
         0.22678912]),
 array([-0.14483408, -1.18666949,  0.52319288, -0.14895357,  1.13640931,
         0.24464938,  0.47165948]),
 array([-0.06464781, -1.14507694,  0.71639249, -0.17427699,  1.22535672,
         0.04290509,  0.44703641, -0.17000883]),
 array([-0.05180161, -1.17885587,  0.7131505 , -0.2668673 ,  1.23852543,
        -0.0112269 ,  0.53356062, -0.16512392,  0.07556192]),
 array([-0.05076509, -1.18112097,  0.72046966, -0.26702131,  1.25551498,
        -0.01488768,  0.54334331, -0.18129495,  0.07485133, -0.01371757]),
 array([-0.04826534, -1.1947611 ,  0.75350698, -0.3660346 ,  1.25822796,
        -0.2436798 ,  0.59200253, -0.31258593,  0.29008666, -0.00446666,
         0.18222971]),
 array([-0.07979969, -1.19398816,  0.70330828, -0.31194246,  1.15578354,
        -0.20151168,  0.37426962, -0.24924464,  0.15969434,  0.20228348,
         0.19058189,  0.17304726]),
 array([-0.07738865, -1.19133281,  0.70612667, -0.30971746,  1.15231086,
        -0.19629704,  0.37146199, -0.23314129,  0.15534809,  0.21208256,
         0.17394624,  0.17193543, -0.01393284]),
 array([-0.0753534 , -1.21644843,  0.68071732, -0.34069761,  1.12961825,
        -0.16224072,  0.31720036, -0.20446703, -0.01297671,  0.25732481,
         0.07079818,  0.3459604 , -0.00262823,  0.14607586]),
 array([-0.05644964, -1.21678855,  0.72548823, -0.33153558,  1.16291879,
        -0.16392005,  0.29074018, -0.16341797, -0.03397237,  0.40350929,
         0.02670832,  0.43405238, -0.16004946,  0.13632434, -0.12941052]),
 array([-0.05827822, -1.21486228,  0.72322674, -0.32540242,  1.16329618,
        -0.15821846,  0.29026015, -0.16572706, -0.02986421,  0.4011931 ,
         0.04314036,  0.42936778, -0.14979831,  0.11913112, -0.13020816,
        -0.01413   ]),
 array([-0.05934892, -1.22472885,  0.73225393, -0.33675343,  1.19583165,
        -0.15494949,  0.32066067, -0.16799003, -0.04242223,  0.42318765,
         0.03115131,  0.5175169 , -0.17445577,  0.17393384, -0.22226471,
        -0.01854605, -0.0757753 ]),
 array([-0.0672958 , -1.22667385,  0.70894407, -0.31851223,  1.17753572,
        -0.10067525,  0.32392765, -0.12360851, -0.04687123,  0.40556981,
         0.06478039,  0.50126667, -0.04904373,  0.13861705, -0.14547006,
        -0.14698867, -0.08199948, -0.10487433]),
 array([-0.05192069, -1.21465231,  0.73049335, -0.29718558,  1.15721375,
        -0.09348519,  0.25043941, -0.13310564, -0.10632983,  0.41244137,
         0.08290202,  0.45377724, -0.03428422, -0.03401566, -0.09877456,
        -0.25092347,  0.09783713, -0.09500843,  0.14660507]),
 array([-0.04741836, -1.21757008,  0.73349799, -0.3048916 ,  1.15418033,
        -0.09452983,  0.24938652, -0.11916986, -0.10378385,  0.4251077 ,
         0.07963657,  0.44968948, -0.02659308, -0.03688665, -0.06323581,
        -0.26005022,  0.12027103, -0.13231115,  0.14501056, -0.03071062]),
 array([-0.0443876 , -1.23188087,  0.7465555 , -0.31676089,  1.17984413,
        -0.08828922,  0.25302678, -0.11654545, -0.14816276,  0.41724854,
         0.03768358,  0.45993169, -0.01483245, -0.06149808, -0.05390686,
        -0.37395385,  0.15036015, -0.20469853,  0.26516999, -0.026031  ,
         0.0986879 ]),
 array([-0.03953258, -1.23316148,  0.75960071, -0.32683117,  1.1872412 ,
        -0.10668613,  0.2503748 , -0.11957089, -0.14889246,  0.43987518,
         0.03953745,  0.48045851, -0.02212142, -0.06723161, -0.04145904,
        -0.3782973 ,  0.20840336, -0.22028179,  0.30189728, -0.08663419,
         0.09650422, -0.04919566]),
 array([-0.0371446 , -1.23784583,  0.76380597, -0.34148539,  1.19793377,
        -0.11680211,  0.26873751, -0.11755845, -0.14562901,  0.44094896,
         0.01621579,  0.47853934, -0.04347314, -0.06000431, -0.03565502,
        -0.39045059,  0.21358194, -0.27791097,  0.3177618 , -0.12350553,
         0.15636239, -0.04727673,  0.04854042]),
 array([-0.02922925, -1.24555512,  0.78930355, -0.36162509,  1.2497503 ,
        -0.16212028,  0.30356572, -0.18122814, -0.15144317,  0.43116423,
         0.00912675,  0.55657342, -0.04082888,  0.01190001, -0.05940233,
        -0.40962052,  0.25740422, -0.29695757,  0.51310552, -0.1791906 ,
         0.2809141 , -0.2491288 ,  0.04248335, -0.16306721]),
 array([-0.04254558, -1.24208586,  0.7689593 , -0.33868521,  1.23511732,
        -0.12021928,  0.2793157 , -0.1602081 , -0.18489343,  0.42631334,
         0.01009852,  0.55323927,  0.00462178,  0.01264532, -0.02419278,
        -0.42198761,  0.24260485, -0.27216791,  0.49986652, -0.07713402,
         0.25138323, -0.18467303, -0.05923064, -0.16545412, -0.08166157]),
 array([-0.0465226 , -1.25014369,  0.76607469, -0.34767902,  1.24736001,
        -0.1239758 ,  0.30365985, -0.17346304, -0.17307825,  0.40576199,
         0.0089203 ,  0.55385511,  0.00484687,  0.03958879, -0.02370097,
        -0.40122559,  0.23360029, -0.27997025,  0.51346956, -0.08298886,
         0.31153506, -0.20116744, -0.02178132, -0.22594532, -0.0837336 ,
        -0.04870131]),
 array([-0.05355462, -1.26223405,  0.73345028, -0.35082405,  1.21831329,
        -0.07899303,  0.29167703, -0.09932279, -0.21350336,  0.43949171,
        -0.04901296,  0.55043291,  0.01056313,  0.04028864,  0.05627058,
        -0.39993758,  0.29218857, -0.30496115,  0.48842311, -0.03914319,
         0.2936341 , -0.02106021, -0.07198295, -0.11533123, -0.26424276,
        -0.05541874, -0.14439073]),
 array([-0.05764967, -1.26380577,  0.72595612, -0.35409494,  1.21627179,
        -0.07959031,  0.30000476, -0.10043292, -0.19965124,  0.43084273,
        -0.04072622,  0.53909032,  0.01215901,  0.04143126,  0.05657016,
        -0.38432681,  0.29079851, -0.29249677,  0.48236796, -0.04196007,
         0.30190632, -0.02330053, -0.03743048, -0.12528092, -0.24344145,
        -0.09121684, -0.14590959, -0.02836091]),
 array([-0.05750069, -1.26303929,  0.7264353 , -0.3528161 ,  1.21692991,
        -0.07939369,  0.30012716, -0.10201889, -0.19943082,  0.42830877,
        -0.03918969,  0.53756271,  0.01417794,  0.04113408,  0.05635251,
        -0.38439068,  0.28796658, -0.29228283,  0.48010467, -0.04091127,
         0.30243391, -0.0248765 , -0.03701238, -0.1316702 , -0.24158133,
        -0.09503041, -0.13927061, -0.02805806,  0.00525317]),
 array([-0.05737803, -1.26369441,  0.72318349, -0.35503495,  1.21128927,
        -0.08246803,  0.29926297, -0.10259973, -0.19236934,  0.42735354,
        -0.02797981,  0.53073825,  0.02090163,  0.03215901,  0.05766827,
        -0.38343025,  0.28829762, -0.27973137,  0.47918964, -0.03091076,
         0.29777744, -0.02725852, -0.03000476, -0.13352395, -0.21316744,
        -0.10326825, -0.1223092 , -0.05754855,  0.00391059, -0.02334883]),
 array([-0.0576237 , -1.26365326,  0.72257798, -0.35632187,  1.2102027 ,
        -0.08471094,  0.29785805, -0.10291543, -0.19265615,  0.43048671,
        -0.02830505,  0.5357802 ,  0.01795834,  0.03519243,  0.05363388,
        -0.38282347,  0.28863599, -0.27951144,  0.48477398, -0.03120516,
         0.30227398, -0.0292826 , -0.0310843 , -0.13037516, -0.21403516,
        -0.09052326, -0.12604481, -0.04993933, -0.00938579, -0.02395255,
        -0.01052184]),
 array([-0.05535703, -1.25849326,  0.72459992, -0.34556364,  1.23735602,
        -0.06520989,  0.34396677, -0.07482924, -0.18595979,  0.43679494,
        -0.0934227 ,  0.5425026 , -0.08647455,  0.09540645, -0.00854579,
        -0.30035337,  0.27708186, -0.28709281,  0.48090528, -0.14662611,
         0.30837162, -0.1220206 ,  0.01041883, -0.1082045 , -0.27820151,
        -0.07227433, -0.38675386,  0.02682164, -0.16504782,  0.24827113,
         0.0018918 ,  0.21542593]),
 array([-0.07660301, -1.25867984,  0.70011464, -0.32928611,  1.23471079,
        -0.027067  ,  0.35109469, -0.04739213, -0.17528832,  0.4357674 ,
        -0.08138864,  0.51209002, -0.07201382,  0.04797805,  0.0197682 ,
        -0.32768006,  0.30670366, -0.28624999,  0.471496  , -0.13809772,
         0.25486831, -0.11280696, -0.03265926, -0.08986456, -0.27082161,
        -0.10619742, -0.38032265, -0.09521031, -0.13096725,  0.1768088 ,
         0.12600837,  0.22088541,  0.09862315]),
 array([-0.07349535, -1.25171964,  0.70408522, -0.32371478,  1.23058395,
        -0.03006712,  0.33911055, -0.05073846, -0.18382203,  0.43293573,
        -0.08241775,  0.50853542, -0.0639828 ,  0.04362653,  0.03462526,
        -0.33669993,  0.31636803, -0.29657534,  0.47211891, -0.13658591,
         0.25259912, -0.09667077, -0.03522385, -0.07613333, -0.27634503,
        -0.10769076, -0.3692595 , -0.0960632 , -0.09206094,  0.16643285,
         0.14806931,  0.18122383,  0.09620935, -0.03151046]),
 array([-0.07220032, -1.25567368,  0.69663722, -0.32980018,  1.22374384,
        -0.02628357,  0.34305859, -0.03556251, -0.17939612,  0.44429305,
        -0.0792888 ,  0.50998306, -0.06000979,  0.03324513,  0.04023871,
        -0.35610323,  0.32855679, -0.30957754,  0.48595671, -0.13800895,
         0.25080614, -0.09404118, -0.05612382, -0.0727461 , -0.29413797,
        -0.10013598, -0.36717423, -0.11000008, -0.09082524,  0.11585789,
         0.16137345,  0.1522871 ,  0.14765295, -0.02848992,  0.04109834]),
 array([-0.06993027, -1.25724731,  0.70479277, -0.32138867,  1.23265723,
        -0.01988421,  0.3380419 , -0.04163832, -0.19967684,  0.43876208,
        -0.09553539,  0.50596496, -0.06310977,  0.0280508 ,  0.05409189,
        -0.3637261 ,  0.3553984 , -0.32667693,  0.50410441, -0.15767816,
         0.25302871, -0.0922049 , -0.05943843, -0.0445774 , -0.29851745,
        -0.07559564, -0.3770831 , -0.11196436, -0.07187654,  0.11440613,
         0.22896643,  0.13407072,  0.18613142, -0.09784654,  0.03711038,
        -0.05523458]),
 array([-7.47534604e-02, -1.25400676e+00,  6.96248627e-01, -3.05135324e-01,
         1.24436454e+00,  1.09569442e-04,  3.48032054e-01, -4.79147146e-02,
        -2.09453775e-01,  4.05834478e-01, -1.02136540e-01,  4.79897859e-01,
        -6.70023499e-02,  2.28605245e-02,  4.60403815e-02, -3.41631153e-01,
         3.41629650e-01, -2.82657586e-01,  4.75578368e-01, -1.26644104e-01,
         2.21267465e-01, -8.74814971e-02, -5.69889835e-02, -5.00882600e-02,
        -2.54335638e-01, -8.39379705e-02, -3.38769571e-01, -1.29400518e-01,
        -7.55124729e-02,  1.43924581e-01,  2.27230101e-01,  2.41708669e-01,
         1.58067156e-01, -3.63027077e-02, -7.26748123e-02, -6.13410270e-02,
        -8.73218782e-02]),
 array([-8.37811232e-02, -1.26034843e+00,  6.88735232e-01, -3.08888434e-01,
         1.26070611e+00,  2.50983168e-02,  3.71523952e-01, -3.30352524e-02,
        -2.17260537e-01,  3.92456568e-01, -1.37159807e-01,  4.71220038e-01,
        -9.32965206e-02,  1.76822126e-02,  4.01486469e-02, -3.50675318e-01,
         3.64505109e-01, -2.95750528e-01,  5.24745441e-01, -1.55866304e-01,
         2.56586417e-01, -1.22800605e-01, -5.22291564e-02, -4.77248534e-02,
        -2.61262592e-01, -3.43243326e-02, -3.49328829e-01, -8.74438296e-02,
        -9.71665885e-02,  1.38970978e-01,  2.63210958e-01,  2.41719997e-01,
         2.86714222e-01, -6.78487403e-02, -6.94021803e-04, -1.90984942e-01,
        -9.50501709e-02, -1.03383745e-01]),
 array([-0.08626541, -1.26263246,  0.68414591, -0.30890511,  1.25907573,
         0.03198798,  0.37733242, -0.02671036, -0.2139211 ,  0.39012168,
        -0.13926106,  0.46282576, -0.09412133,  0.01140414,  0.03900183,
        -0.35193037,  0.36155424, -0.28958482,  0.52100001, -0.14325681,
         0.24947961, -0.11404164, -0.06065579, -0.04676009, -0.26083769,
        -0.03656622, -0.33800554, -0.09073974, -0.08773596,  0.13375027,
         0.26241713,  0.25064762,  0.28731733, -0.03755431, -0.00811653,
        -0.17443482, -0.12533601, -0.10539698, -0.02402974]),
 array([-0.08744464, -1.26780471,  0.67799518, -0.31746532,  1.25867741,
         0.03014504,  0.39143222, -0.01441009, -0.20104325,  0.39668534,
        -0.14356661,  0.4583728 , -0.1107086 ,  0.00960969,  0.02620149,
        -0.35422507,  0.35857762, -0.2951813 ,  0.53324297, -0.15028699,
         0.27504715, -0.12825272, -0.04291289, -0.06403071, -0.25892372,
        -0.03600658, -0.34262444, -0.06802704, -0.09457005,  0.15289509,
         0.25191917,  0.24933684,  0.30583453, -0.03598453,  0.05367132,
        -0.18959402, -0.09176225, -0.16735937, -0.02826312, -0.04907397])]
In [56]:
#checking whether the coefficient are correct
#should preserve the acf!
for p in range(1, len(coef)+1):
    computeACF(signalAfterHamming, 1) == coef[0][0]*computeACF(signalAfterHamming, 0)
    Rind = np.array([[ np.abs(i-k) for k in range(1, p+1)] for i in range(1, p+1)])
    R = np.array([[computeACF(signalAfterHamming, np.abs(i-k)) for k in range(1, p+1)] for i in range(1, p+1)])
    acfArr = np.array([computeACF(signalAfterHamming, i) for i in range(1, p+1)])
    coefArr = coef[p-1]
    if np.sum((np.round(np.matmul(R, coefArr) - acfArr, 7))**2)==0:
        print(f"Order {p} coefficients verified and correct...proceed")
    else:
        print("Error <Coefficients don't match>")
        break
Order 1 coefficients verified and correct...proceed
Order 2 coefficients verified and correct...proceed
Order 3 coefficients verified and correct...proceed
Order 4 coefficients verified and correct...proceed
Order 5 coefficients verified and correct...proceed
Order 6 coefficients verified and correct...proceed
Order 7 coefficients verified and correct...proceed
Order 8 coefficients verified and correct...proceed
Order 9 coefficients verified and correct...proceed
Order 10 coefficients verified and correct...proceed
Order 11 coefficients verified and correct...proceed
Order 12 coefficients verified and correct...proceed
Order 13 coefficients verified and correct...proceed
Order 14 coefficients verified and correct...proceed
Order 15 coefficients verified and correct...proceed
Order 16 coefficients verified and correct...proceed
Order 17 coefficients verified and correct...proceed
Order 18 coefficients verified and correct...proceed
Order 19 coefficients verified and correct...proceed
Order 20 coefficients verified and correct...proceed
Order 21 coefficients verified and correct...proceed
Order 22 coefficients verified and correct...proceed
Order 23 coefficients verified and correct...proceed
Order 24 coefficients verified and correct...proceed
Order 25 coefficients verified and correct...proceed
Order 26 coefficients verified and correct...proceed
Order 27 coefficients verified and correct...proceed
Order 28 coefficients verified and correct...proceed
Order 29 coefficients verified and correct...proceed
Order 30 coefficients verified and correct...proceed
Order 31 coefficients verified and correct...proceed
Order 32 coefficients verified and correct...proceed
Order 33 coefficients verified and correct...proceed
Order 34 coefficients verified and correct...proceed
Order 35 coefficients verified and correct...proceed
Order 36 coefficients verified and correct...proceed
Order 37 coefficients verified and correct...proceed
Order 38 coefficients verified and correct...proceed
Order 39 coefficients verified and correct...proceed
Order 40 coefficients verified and correct...proceed
In [57]:
#computing the minima energy directly from formula
Emin = []
for p in range(1,len(coef)+1):
    acfArr = np.array([computeACF(signalAfterHamming, i) for i in range(1, p+1)])
    E = computeACF(signalAfterHamming, 0) - np.sum(acfArr*coef[p-1])
    Emin.append(E)
In [58]:
error = []
for p in range(1, len(coef)+1):
    signalShifted = np.pad(np.array(signalAfterHamming, dtype = float),(p, 3*len(signalAfterHamming)-p), 'constant')
    signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(0, 3*len(signalAfterHamming)), 'constant')
    signalEstimate = np.zeros_like(signalShifted)
    for i in range(0, len(signalAfterHamming) + p-1):
        signalEstimate[i] = np.sum(signalShifted[i:(i+p)]*coef[p-1][::-1])
    error.append(np.sum((signalToEstimate - signalEstimate)**2))
In [59]:
plt.figure(figsize=(16, 8))
plt.gcf().set_dpi(300)
plt.suptitle("Error computed for different orders of LPA for different methods of error computation")
plt.subplot(121)
plt.plot([i for i in range(1, order + 1)], error, 'ro')
plt.plot([i for i in range(1, order + 1)],  error, 'b', alpha = 0.3, linewidth = 1)
plt.vlines([i for i in range(1, order + 1)] , 0, error, color='g', linestyles='--', alpha = 0.3)
plt.ylim(bottom=0)
plt.title("Calculated MSE error/loss from the computed coefficients")
plt.xlabel("Order of LPA (p)")
plt.ylabel("Calculated error")
plt.legend(["Error computed"])
plt.subplot(122)
plt.plot([i for i in range(1, order + 1)], Emin, 'ro')
plt.plot([i for i in range(1, order + 1)], Emin, 'b', alpha = 0.3, linewidth = 1)
plt.vlines([i for i in range(1, order + 1)] , 0, Emin, color='g', linestyles='--', alpha = 0.3)
plt.ylim(bottom=0)
plt.title("Calculated MSE error/loss from direct loss equation")
plt.xlabel("Order of LPA (p)")
plt.ylabel("Cmputed error from formula")
plt.legend(["Error computed from formula"])
Out[59]:
<matplotlib.legend.Legend at 0x7fd580ea1e90>

Order of lpa p=6

In [60]:
z, p, k =tf2zpk(np.sqrt(Emin[5]), np.insert((-1)*coef[5] , 0, 1))
In [61]:
#Zeros
z
Out[61]:
array([], dtype=float64)
In [62]:
#poles
p
Out[62]:
array([ 0.96823923+0.j        ,  0.10438044+0.97037854j,
        0.10438044-0.97037854j, -0.46153749+0.7935437j ,
       -0.46153749-0.7935437j , -0.29179197+0.j        ])
In [63]:
angle = np.linspace( 0 , 2 * np.pi , 150 ) 

radius = 1

x = radius * np.cos( angle ) 
y = radius * np.sin( angle )
plt.figure(figsize=(10,10)) 
plt.plot(x, y, 'r', alpha = 0.2)
plt.plot([i.real for i in p], [j.imag for j in p], 'x', label = "Poles")
plt.xlim([-1.1, 1.1])
plt.ylim([-1.1, 1.1])
plt.title("Pole zero plot for LPA of order 6")
plt.xlabel("Real axis")
plt.ylabel("Imaginary axis")
plt.legend()
Out[63]:
<matplotlib.legend.Legend at 0x7fd580dd10d0>

Order of lpa p=10

In [64]:
z, p, k =tf2zpk(np.sqrt(Emin[9]),np.insert((-1)*coef[9] , 0, 1))
In [65]:
#Zeros
z
Out[65]:
array([], dtype=float64)
In [66]:
#poles
p
Out[66]:
array([ 0.98596008+0.j        , -0.60881571+0.67251047j,
       -0.60881571-0.67251047j,  0.12965731+0.9607452j ,
        0.12965731-0.9607452j , -0.24424742+0.81540094j,
       -0.24424742-0.81540094j,  0.09709139+0.32491389j,
        0.09709139-0.32491389j,  0.21590368+0.j        ])
In [67]:
angle = np.linspace( 0 , 2 * np.pi , 150 ) 

radius = 1

x = radius * np.cos( angle ) 
y = radius * np.sin( angle )
plt.figure(figsize=(10,10)) 
plt.plot(x, y, 'r', alpha = 0.2)
plt.plot([i.real for i in p], [j.imag for j in p], 'x', label = "Poles")
plt.xlim([-1.1, 1.1])
plt.ylim([-1.1, 1.1])
plt.title("Pole zero plot for LPA of order 10")
plt.xlabel("Real axis")
plt.ylabel("Imaginary axis")
plt.legend()
Out[67]:
<matplotlib.legend.Legend at 0x7fd580d50cd0>
In [68]:
#can we obtained by sqrt(Emin)

for i in range(0, order):
    print(f"Gain for order {i+1} lpa is {np.sqrt(Emin[i])}")
Gain for order 1 lpa is 0.36339828510963124
Gain for order 2 lpa is 0.24873498775341477
Gain for order 3 lpa is 0.2403009620959485
Gain for order 4 lpa is 0.2391857806184903
Gain for order 5 lpa is 0.15177843115736972
Gain for order 6 lpa is 0.14782367504819888
Gain for order 7 lpa is 0.13034807468633497
Gain for order 8 lpa is 0.12845053769737838
Gain for order 9 lpa is 0.1280833119118877
Gain for order 10 lpa is 0.1280712605166856
Gain for order 11 lpa is 0.12592683354762932
Gain for order 12 lpa is 0.12402704108715562
Gain for order 13 lpa is 0.12401500218105116
Gain for order 14 lpa is 0.12268474172919244
Gain for order 15 lpa is 0.12165309845022389
Gain for order 16 lpa is 0.12164095340764093
Gain for order 17 lpa is 0.121291225807064
Gain for order 18 lpa is 0.12062236316085043
Gain for order 19 lpa is 0.11931905098427194
Gain for order 20 lpa is 0.11926277030476519
Gain for order 21 lpa is 0.11868058129175943
Gain for order 22 lpa is 0.11853687815391492
Gain for order 23 lpa is 0.11839714915354366
Gain for order 24 lpa is 0.11681239888781607
Gain for order 25 lpa is 0.11642225905285154
Gain for order 26 lpa is 0.11628411093351615
Gain for order 27 lpa is 0.11506554073484086
Gain for order 28 lpa is 0.11501925546157239
Gain for order 29 lpa is 0.11501766842778051
Gain for order 30 lpa is 0.11498631219339239
Gain for order 31 lpa is 0.11497994700508478
Gain for order 32 lpa is 0.11228023897444306
Gain for order 33 lpa is 0.11173285638681872
Gain for order 34 lpa is 0.11167737232920139
Gain for order 35 lpa is 0.11158301682580743
Gain for order 36 lpa is 0.11141267476192156
Gain for order 37 lpa is 0.11098709496332554
Gain for order 38 lpa is 0.11039237551089325
Gain for order 39 lpa is 0.11036049907255603
Gain for order 40 lpa is 0.11022753086462132

LPC spectra magnitude plots

For order p = 2

In [69]:
p = 2

signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
    signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])

w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]),  0, 1), None, 1)

plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))  

plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[69]:
<matplotlib.legend.Legend at 0x7fd580c7d510>

For order p = 4

In [70]:
p = 4

signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
    signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])

w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]),  0, 1), None, 1)

plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))  

plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[70]:
<matplotlib.legend.Legend at 0x7fd580c07c90>

For order p = 6

In [71]:
p = 6

signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
    signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])

w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]),  0, 1), None, 1)

plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))  

plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[71]:
<matplotlib.legend.Legend at 0x7fd581111990>

For order p = 8

In [72]:
p = 8

signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
    signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])

w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]),  0, 1), None, 1)

plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))  

plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[72]:
<matplotlib.legend.Legend at 0x7fd580bef690>

For order p = 10

In [73]:
p = 10

signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
    signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])

w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]),  0, 1), None, 1)

plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))  

plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[73]:
<matplotlib.legend.Legend at 0x7fd580aeeb90>

While reconstruction, we would now model the error as white noise since the sound is unvoiced we no longer have the glottal vibrations

In [74]:
signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(0, 10*len(signalAfterHamming)), 'constant')
signalReconstructed = np.zeros_like(signalToEstimate)
mean = 0
std = 1 
num_samples = len(signalToEstimate)
whiteNoise = np.random.normal(mean, std, size=num_samples)
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(whiteNoise)
#plt.vlines(np.arange(0, signalToEstimate.shape[0], 59), 0, 1, colors = 'r', linestyles='--')
plt.xlabel("Times(samples)")
plt.ylabel("Amplitude")
plt.title("White noise generated")
Out[74]:
Text(0.5, 1.0, 'White noise generated')
In [75]:
P = 10#choosing the order of analysis
G = np.sqrt(Emin[P-1])
whiteNoise = np.pad(np.array(whiteNoise, dtype = float),(P, P), 'constant', constant_values=0)
#signal y[n] is assumed casual
for i in range(P, signalToEstimate.shape[0]-2000):
    signalReconstructed[i] = G*whiteNoise[i] + np.sum(signalReconstructed[i-P:i]*coef[P-1][::-1])
In [76]:
plt.figure(figsize=(16,8))
plt.gcf().set_dpi(300)
plt.suptitle("Signal estimated along with ground truth")
plt.subplot(121)
plt.plot([i*(1.0/16000) for i in range(len(signalReconstructed))][0:480],signalReconstructed[0:480])
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Estimate of the signal(after resynthesis and de-emphasis)")
plt.subplot(122)
plt.plot([i*(1.0/16000) for i in range(len(signalReconstructed))][0:480],signalToAnalyse[0:480])
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Ground truth(true signal")
Out[76]:
Text(0.5, 1.0, 'Ground truth(true signal')
In [77]:
scaled = np.int16(signalReconstructed/np.max(np.abs(signalReconstructed)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('ssReconstructed.wav', 16000, scaled)
In [77]:
 
In [77]: